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1.  INTRODUCTION 


Accurate  knowledge  of  electron  densities  in  die  ionosphere  is  paramount  for  enabling  precise  navigation, 
uninterrupted  ground-to-satellite  communicatians,  high  accuracy  positioning,  and  surveillance.  Thus,  an  ability  to 
forecast  regional  ionospheric  conditions  has  applications  in  most  areas  of  military  operatirms,  fiom  HF 
communications  to  delivering  vital  data  to  individual  soldiers  in  remote  battlefields,  to  predsion-guided  weapons, 
to  space  based  intelligence  gathering  and  space  object  tracking. 

During  the  first  year  of  the  AF  sponsored  SBIR  Phase  It  investigation  we  developed  a  new  computer  model  for 
simulating  time  evolution  of  ion  and  electron  densities  in  the  iono^here  on  a  global  scale  and  developed  the 
corresponding  tangent  linear  and  adjoint  models  with  respect  to  parallel  transport 

These  new  computational  tools  will  further  evolve  into  a  prototype  numerical  ionospheric  forecast  systems,  that 
will  assimilate  available  measurements  of  ionospheric  electron  cont^t  firom  networks  of  ground-based  GPS 
reference  stations  and  otiier  instruments. 

Significant  improvements  in  the  reliability  of  convaitional  weatibier  forecasts  in  the  last  decade  or  so  are  largely 
due  to  advances  in  data  assimilation  techniques.  The  atmosphere,  including  the  ionosphere,  is  a  chaotic  system; 
small  errors  in  the  mitial  conditions  of  a  forecast  grow  rapidly,  and  affect  predictability.  Furthermore,  predictability 
is  limited  by  model  errors  due  to  the  approximate  simulation  of  relevant  physical  processes  in  the  numerical 
models  and  to  a  poorly  known  external  forcing.  Data  assimilation  aims  to  decrease  these  uncertainties  by  using 
observations  to  obtain  better  initial  conditions  and/or  to  provide  better  estimates  of  poorly  known  err^)irical 
quantities  in  parameteiizations  of  various  physical  processes  in  the  models. 

Development  of  practical  and  effidoit  data  assimilation  schemes  as  well  as  operational  implementations  of 
forecast  systems  depend  aitically  on  the  design,  quality,  and  maintainability  of  the  underlying  physical  propagator 
model. 

In  the  course  of  the  Phase  I  investigation  it  was  decided  that  the  computer  code  of  the  model  available  to  us. 
Coupled  Thermosphere  Ionosphere  Model  (ClIM),  is  unacceptable  for  the  purpose  of  adding  on  data  assimilation 
capabilities  and  practical  implementation  of  the  forecast  system  for  several  reasons,  with  major  reasons  listed 
below: 

-  Very  poor  quality  of  code  leading  to  inefficiencies,  redundancy,  and  making  code  hard  to  understand  and 
modify, 

-  Almost  complete  lack  of  documentation. 

-  Extensive  use  of  Fortran  77  and  Fortran  66  features  that  are  being  phased  out  in  foe  modem  and  fixture  compiler 
versions. 

The  first  major  task  of  the  Phase  n  investigation  was  foen  creating  a  new  computer  code  implementing  foe  model 
and  performing  foe  first  round  of  validation.  This  task  is  largely  completed,  certain  features  of  the  model  that  are 
deemed  as  non-aitical  are  stiU  left  in  foe  developmait  stage  (.e,g.,  correct  calculations  of  ni^t  time  ionization 
rates;  high-latitude  transport)  and  will  be  added  afiier  testing  foe  model  and  assimilation  scheme  with  real  GPS 
data. 

While  development  of  foe  propagator  model  “from  scratch”  took  time,  we  deem  it  to  be  an  absolutely  cmcial 
component  that  will  ^able  and  strongly  facilitate  reaching  the  end  objective.  For  example,  building  foe  tangent 
linear  and  adjoint  models  foat  are  necessary  for  implementing  data  assimilation  capabilities  for  existing  numerical 
global  atmospheric  models  usuaEy  take  several  person-years.  In  this  Phase  H  effort  foe  development  of  bofo 
components  vrifo  respect  to  parallel  transport  took  less  foan  2  person-days. 
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We  believe  that  these  reductions  in  the  development  time  are  due  to  adherence  to  modem  software  development 
practices  and  tools,  object-oriented  programming  style,  and  maintaining  well-documented  code  throughout  the 
durati(m  of  the  project 

The  Mowing  personnel  contributed  to  the  development  of  the  model:  Dr.  Boris  Khattatov,  PI;  Dr.  Michael 
Murphy,  Semor  Computer  Scientist;  Mr.  James  Adams,  Senior  Software  Engineer,  Dr  Timothy  Fuller-RowelL, 
Consultant 


2.  MODEL 

The  developed  model  is  a  numerical  global  model  of  the  ionosphere  system  loosely  based  on  Millward  et  al 
(1996),  Bailey  and  Balan  (1996),  FuUer-RoweU  (1996)  and  Huba  et  al  (2000). 

The  dynamic  equations  and  vertical  ExB  transport  for  seven  ions  (H,  O,  O2,  He,  NO,  N2,  N)  are  solved  on  a  fixed 
Euletian  grid  in  magnetic  p,  q,  and  longitude  coordinates.  An  exatrq^le  of  the  low-latitude  configuration  is  shown 
below  (only  20  bngitudes  and  30  p  values  are  shown  for  clarity,  regular  model  configuration  is  100x100x100). 
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Figure  1.  A  part  of  the  model  magnetic  grid. 


The  model  solves  plasma  dynamics  equations  -  parallel  and  ExB  continuity  and  momentum  —  for  seven  ion 
species  and  electrons  and  energy  conservation  equation  for  the  three  major  ions  and  electrons.  The  model  includes 
chemical  interactions  with  neutrals  and  ion-ion  and  ion-neutral  collision  rates  and  Photoionization.  ExB  drift  is 
computed  by  via  Fejer&Scherless  model  at  the  equator  and  by  Weimer  (2000)  model  at  high  latitude.  The  high- 
latitude  transport  is  currently  tumed-off  pending  model  validation  at  low  latitudes. 

Model  prognostic,  external  and  diagnostic  variables  are  listed  in  Tables  1-3. 
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Variable 

Units 

JH^ription 

Comine^ 

^atne 

Ion  densities 

particles/ 

Local  (pohn)  volume  density  of  a  particular 

at  psresent  there  are  7 

ion  ^secies 

ions:  o',  H',  He',  N,', 
Oj‘,NO',N' 

Ion 

tempexatuxes 

Lo(ail  tempers 
if^jecies 

ntne 

Ion  velocities 

m/s 

Local  ^hxt)  velocity  of  a  pslictilar  ion 
species  akmg  Ihe  magnetic  field  line 
passing  through  tiiis  point 

same 

Electron 

tisn|)er8tuie 

K,  degree 

Loed  deeta  ten^peratiae 

Hectron 

vdocity 

m/s 

Local  (point)  velocity  of  electrons  along  the 
magnetic  field  line  passing  tiirough 
point 

Blectixm 

partictes/ 

Local  (^hrt)  viilisme  tUmsity  of  electro 

^  sum  <^atl  local 

ns* 

iondetisiti^ 

Table  1.  Model  Prognostic  Variables 


Varbibte 

Nantt 

Units 

Pmiiptioii 

Omnwnts 

Neutral 

particles/ 

Local  (point)  volume  density 

-  at  present  there  axe  7  neutrals: 

densities 

m’ 

of  a  particular  neutral  impedes 

0,0z,N2,He,H,N0,N. 

Neidnil 

K. 

Loed  tempetatine 

aanie 

tempecaiise 

deftee 

iieidxa}  i^des  (one  fin  ali) 

Neutral  20(Qal 

m/s 

Local  (point)  velocity  of  all 

same 

vdocity 

neutral  spedes(Qne  fin  a&) 
in  tile  zonal  direction  (east- 
west,  eastward  is  positive) 

Neutral 

Local  (point)  vdocity  of  d! 

inen^istal 

neutral  ^pec^  (tme  for  all} 

vdoetty 

in  tiie  meddiond  direction 
(nortihsouih^  nottiiward  is  ; 
positive) 

Table  2.  Model  External  Variables 
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Variable  Name 

Units 

Dcicripflon 

Comments 

ExB  zonal  vclodty 

m/s 

Local  (point)  velocity  assodated 
with  the  zonal  ExB  drift  of  the 
magnetic  field  line  passing  through 
this  point 

Needs  to  be  computed  from  empirical  ExB 
models 

ExB  t&eridional 
vdocity 

ffi/S 

Local  (point)  velocity  assodated 
with  the  meridional  ExB  drift  of  the 
magnetic  field  line  passing  throu^ 
this  point 

same  as  above 

Photo  production 

Partides/s/ra^ 

Number  of  particles  of  a  particular 
ion  species  produced  as  a  result  of 
photoionization  per  second  per  unit 
volume. 

-  at  present  there  are  7  neutrals,  only  5  of  those 
can  be  photoionized:  0, 02,  Nj,  He,  N.  Other 
ions  are  produced  via  chemical  reactions,  such 
as 

O'+H  -♦H'  +0. 

Ch<^lcal 

production 

Partides/s/tc* 

Number  of  particles  of  a  particular 
icm  spedes  produced  as  a  result  of 
chemical  reactions  per  second  per 
unit  volume. 

There  axe  21  chemical  reactions  at  the  present 
e.g.,0'+H  -*H'  +0. 

Chemical  loss 

Partides/s/m^ 

Number  of  particles  of  a  particular 
ion  species  lost  as  a  result  of 
chemical  and  recombination 
reactions  per  second  per  unit  volume 

This  value  is  a  product  of  the  density 
(concentration)  of  the  ion  species  being 
destiucted  and  the  chemical  hsx  rate,  L. 

Photoionization 

rates 

l/s 

Coefficients  needed  to  compute 
photo  production 

-  at  present  there  are  7  neutrals,  only  5  of  those 
can  be  photo-iemized:  0,  Q2,  Nj,  He,  N. 

Chemical  reaction 

rates 

mVs 

Coefficients  needed  to  compute 
chemical  loss  due  to  electron 
exchange  reactions. 

There  are  21  chemical  reactions  at  the  present 
c.g.,0'  +H  ->H'  +0. 

Recombination 
reaction  rates 

1/s 

Coeffidents  needed  to  compute  loss 
due  to  recombination  chemical 
reactions 

There  are  7  recombination  reactions,  c.g.,  0’  + 
e  -+0. 

e  represents  an  dectroa 

lon^eutxal 

collision 

frequencies 

1/s 

Drag  on  a  particular  ion  particle  due 
to  collisions  with  a  neutral  i^>ecies. 

There  are  7  ions  and  7  neutrals,  therefore  h  is 
a  7x7  matrix  with  zero  diagonal. 

lon^ion  collision 
froquendes 

]/s 

Drag  on  a  particular  ion  particle  due 
to  collisions  with  a  difterent  ion 
species. 

There  are  7  ions,  therefore  h  is  a  7x7  notrxx 
with  zero  diagonal. 

Ion  heating  rates 

Jhon  thermal 
conductivities 
Electzxm  heating 
rates 

Electron  thermal 
conductivities 

J/mVs 

i/K/m/s 

J/mVs 

J/KAd/s 

Heating  due  to  Joule  heating, 
frictional  collisions  and  other 
processes. 

Heating  due  to  Joule  heating, 
frictional  collisions  and  other 
processes. 

Is  only  computed  fi>r  three  m^or  ions,  O*,  H\ 
He' 

Is  cmlyconqmted  for  three  xsajm  ions,  0\H', 
He‘ 

Table  3.  Model  Diagnostic  Variables 


4 


2.1  Prognostic  Equations 


A  prognostic  equation  allows  one  to  estinuite  a  particular  prognostic  variable  at  a  future  time.  Ibe  prognostic 
variables  are  density,  velocity  and  tempaature  for  ions  and  electron  density,  temperature  and  velocity. 

These  equations  are  given  in  dipole  coordinates,  along  magnetic  flow  tubes.  Therefore  there  is  only  one  dependent 
spatial  coordinate  corresponding  to  the  position  along  the  magnetic  flow  tube.  This  can  be  a  non-chiiiensioiial 
variable  q  or  a  dimensiQnal  variables  =  (Re  is  tbe  radius  of  the  earth). 


211  ConfmoityEqualKmForEadilimS^pedes 

Numerical  solution  of  this  equation  stK)uld  generate  ion  density  (t  4^  t)  given  all  related  variable  at  time  t. 


dN, 

dt 


NV, 


ds 


-fiV,  •  VKl +ViV, -A, 


(1) 


where 


Nf  -density  of imii 

Vj  -velocity  (aligned  with  the  magnetic  flow  tube)  of  ion  i 
5  = 


-diemical  production + photochemical  production 
Lj  -  chemical  loss  rate 
L.  •  -  chetnical  loss 


The  term 


VF, 


6»F^  sin^{eccLat)(}^cos^{eccLat)) 
/?  •  '  (1  +  3  ■  cos^  {eccLat)f 


(2) 


is  a  div^gence  of  ExB  velocity  in  the  vertical  (and  mendional)  plane,  i.e.,  in  p  direction.  is  the  value  of  ExB 

meridional  drifl.  at  the  magnetic  equator  oorre^3ondmg  to  a  particular  p. 


212  MomentinnEqiiatkmForEadifoii^iecies 

Numerical  solution  of  this  equation  should  generate  ion  velodty  given  all  related  variables  at  time  t. 


V,=^ 


i  N  _  Neutrals  N _Ions 

E  £  ''u 

/=1 


ds 


yNi  ds  N,  ds 


•  T  . 

-gsm/  +  -^ 

OT, 

iV  _  Neutrals  . 

+  S  Vj„(F;,cosD-l7„sinD)cos/+  Yj 

a=l  /=*l 


N  _  Ions 


(3) 
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v^ere 


Nf  -density  of  ion  i 

-  velocity  (aligned  with  the  magnetic  flow  tube)  of  ion  i 
Vj  -  velocity  of  ion  j. 

s  =  qK, 

6,  =  -y/l + 3  cos^  {eccLat)  •  | - - 

\eccRadius 

m.  -  mass  of  ion  i. 
k  -  Boltzmann’s  constant. 

7]  -temperature  of  ion  L 

7^  -  electron  temperature 
-  electron  density. 

-zonal neutral  velocity. 

-meridional neutral  velocity. 


-  ion-neutral  collision  frequency. 

-  ion-ion  collision  frequency. 


g  -  acceleration  of  gravity. 

I — inclination  angle  for  this  flow  tube: 
2cos(eccLaO 

yj\+3  CQS^{eccLat) 
sm(eccLat) 
-i-3cos^(eccLat) 

D  -  declination  angle  for  this  flow  tube. 


2J3  EoagyEkfuationForEadilonSpedes 

Numerical  solution  of  tiiis  equation  should  generate  ion  tenperature  givai  all  related  variables  at  time 

t. 


2  ‘I,  a 


+  F,V7:. 


j 


^kN,Tp] 


d(V^ 


ds 


-kNiT.^VVj^+b^—^ 

ds 


dT\  3 


dT 


(4) 

os  J  2  os 


WhCTe 


QandF-  are  the  heating  rates 
^  -  is  the  thermal  conductivity 
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21.4  £lectr(mTat]p^*ati]re  Equation 


It  is  Rtmtlar  to  ion  temperature  equation,  except  that  the  conductivities  and  heating  rates  are  confuted  for 
electrons. 

2.15  Electnm  Density  Equation 

NundferOflnns 

N.=  Z 

/=1 


2JL.6  Electron  Velocity  Equation 
Assumes  tiiat  tiiere  are  no  field-aligned  currents: 


NumberOfU)ns 

I  KK 


(6) 


2.2  Timing 

On  a  Pentium  Xeon  2.2GHz  computer,  1  model  time  step  requires  35  seconds  in  a  configuration  with  10^  grid 
points  and  tire  memory  footprint  is  L5  Gb. 


3.  SELECTED  RESULTS 

In  this  section  we  will  present  selected  results  obtained  from  different  numerical  experiments. 
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Figure  4.  An  example  of  3-D  instantaneous  iso-surface  of  constant  electron  density.  The 
radial  direction  is  exaggerated  to  show  detail. 


Figure  5.  The  x,  y,  and  z-components  of  the  Earth’s  magnetic  field  at  the  North  Pole 
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Figure  6.  The  x-  and  y-component  of  the  Weimer  electric  potential  at  the  North  Foie. 


Figure  7.  Calculated  x-,  y-,  and  z-  component  of  the  ExB  drift  velocities  near  the  North  Pole. 


Figure  8.  Trajectories  of  particles  in  the  magnetic  field  acting  under  ExB  drift  with  different 
starting  locations  near  the  Noth  Pole. 
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4.  TANGENT  LINEAR  AND  ADJOINT  MODELS 


In  Ifae  case  of  the  ionosphere,  lets  assume  that  the  model  M  is  a  forward  time-dependent  discrete  propagator  that 
accepts  the  current  state  i*  (electron  or  ion  densities  throughout  the  model  domain  arranged  in  a  vector)  and  values 
of  several  atmospheric  drivers  pt  (e.g.,  level  of  solar  activity,  etc)  as  inputs  and  generates  state  estimates  for  a  later 
time: 


(7) 

Both  quantities,  Xtand  pt,  are  considered  to  be  model  parameters. 

Generally  speaking,  the  tangmt  linear  of  the  model  M  is  simply  a  derivative  of  the  results  with  respect  to  the  initial 
conditions  or  input  parameters.  Note  tiiat  M  is  a  non-linear  vector  function  and  therefore  its  linearizations  (first 
derivatives)  are  matrices: 


^  dM  dM 

L^  = - and  L„  = - 

"  5x  ^  dp 


(8) 


The  adjoint  of  the  propagator  model  is  simply  a  transposed  of  the  matrix  L  or,  a  way  to  compute  a  product  of  the 
transposed  and  an  arbitraiy  vector.  The  linearization  matrix  describes  sensitivity  of  tiie  model  with  respect  to  the 
initial  conditions  and  the  adjoint  is  used  to  either  minimize  the  cost  function  in  the  variational  assimilation 
approach  or  to  solve  the  Kalman  filter  equations.  Both  these  matrices  are  fundamental  to  implementing  data 
assimilation  schemes. 


These  linearization  matrices  can  be  obtained  in  several  ways: 

-  approximated  via  finite-diffaences  calculations,  Le.,  by  introducing  small  changes  in  Xt  and/or  p*  and  computing 
resulting  changes  in 


^  AM  _ 

AX/  Ax^ 


(9) 


-  by  differentiating  the  actual  computer  code  (e.g.,  Fortran  or  C)  implementing  the  model  and  generating  computer 
code  for  direct  computation  of  L. 

-  by  analytical  differentiation  of  theoretical  equations  of  the  model  and  coding  the  results. 

The  first  approach  can  be  extremely  CPU  intensive  but  is  straightforward  to  implement.  The  other  two  are  much 
more  efficient  but  can  be  bard  to  implement  and  will  have  to  be  re-done  if  changes  are  introduced  into  the  model. 

We  followed  a  variation  of  the  second  approach  and  obtained  and  coded  explicit  calculations  of  both  these 
matrices  for  each  plasma  tube  in  the  model.  An  example  of  the  linearization  matrix  is  shown  in  Figure  9  for 
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Figure  9.  An  example  of  logarithm  of  absolute  value  of  linearization  matrix  for  one  plasma 
tube. 

Ihe  particular  plasma  tube  has  100  points  and  changes  in  ion  densities  at  the  next  time  step  can  be  obtained  by 
inultiplyinglhk  matrix  by  a  vector  of  changes  in  ion  density  for  this  tube  at  aprevious  tie  step.  The  asymmetries  in 
the  matrix  reflect  the  flict  that  in  this  particular  case  tiie  prevailing  parallel  transport  is  directed  from  the  top-left 
end  of  the  tube  to  the  hottom-ri^t 

Once  the  linearization  matrix  is  obtained,  evolution  of  tiie  error-covariance  matrices  for  different  points  on  the 
same  tube  can  be  trivially  computed  as  follows  (IQiattatov  et  al.,  1999): 


C(f+A/)  =  L,C(0^ 


5.  CONCLUSIONS 


The  initial  Ibase  n  development  effort  led  to  oreation  of  a  new  ionospheric  numerical  ^obal  model  by  its  design 
particularly  well  suited  ftir  use  with  a  data  assimilation  scheme.  The  correspondmg  linear  and  adjoint  model  with 
respect  to  parallel  transport  were  iirplemented  and  are  ready  for  exp^iments  with  real  data,  hi  the  remaining  12 
months  of  this  investigation  we  plan  to  implement  a  prototype  system  for  now-casting  and  forecasting  three- 
dimensional  global  electron  doisities  in  tiie  ionosphao  constrained  by  a  sequential  data  assimilation  scheme  and 
validate  the  system  performance. 
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